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ABSTRACT 



Many meteorological forecast applications require the use of grids 
that have a high resolution in a particular area of interest ^ while 
allowing coarser resolution elsewhere. Conventional finite difference 
models often use nested grids to this end. In recent years, finite 
element models have been offered as an alternative. In this study, the 
two-dimensional advection equation with diffusion is defined over a 
rectangular domain. The Galerkin technique is applied to linear basis 
functions on triangular elements. The model is tested to determine the 
sensitivity of the forecast to various nodal geometries. Both equi- 
lateral and right triangular elements are tested. It is found that the 
equilateral arrangement consistently yields a superior forecast. Other 
tests are conducted in which the resolution is varied smoothly versus 
abruptly over the domain. The smoothly varying case gives results that 
are dramatically improved over the abruptly varying case. Among the 
conclusions is the fact that, for a given maximum resolution, the more 
slowly and smoothly the element size is changed, the better the forecast 
obtained. 
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I. INTRODUCTION 



Sub-Synoptic scale meteorological elements are usually not repre- 
sented or forecast well by coarse mesh numerical models. This situation 
may be improved by increasing the resolution of the grid. However^ a 
uniform reduction of the grid mesh requires a significant increase in 
computational time and computer storage. Hence this is a limitation on 
the practical resolution, and therefore accuracy, for a uniform grid. 

The conventional solution is to nest grids so that the geographical area 
or meteorological feature of interest is covered by a fine mesh, and 
less interesting areas or features are left with a coarse mesh* The 
advantage of this technique is to increase accuracy in the desired area 
while keeping the computational storage and time requirements within 
reasonable limits. There are, however, two disadvantages to nested 
grids. First, the abrupt change in grid size at the boundaries of the 
fine mesh generates noise in the solution. Second, boundary conditions 
for the fine mesh must be interpolated from the coarse mesh. 

A better technique would allow the grid to vary smoothly and con- 
tinuously over the domain rather than in abrupt jumps. However, such a 
highly variable grid greatly complicates a finite difference numerical 
model. As an alternative to such a finite difference technique, consider 
the Finite Element Method (FEM) . This method has long been used in 
mechanical engineering but has been adapted to meteorology only during 



11 



the last decade. Pioneering work in meteorological applications was 
done by Cullen (1973) . Staniforth and Mitchell (1978) and Staniforth 
and Daley (1977) have devised more recent forecast models. The most 
recent FEM meteorological model at the Naval Postgraduate School was 
written by Kelly and Williams (1976) and is the precursor to this study. 

The great attractions of the FEM are the more accurate phase speeds 
[Haltiner and Williams (1980)], and its suitability to non-uniform grids. 
The method involves the division of the domain into a number of elements 
and then defining an equation for each node of each element. Thus, 
solution requires solving a large number of linear equations simul- 
taneously. However, there is no limitation on the size of the elements 
or the shape of the domain. The method is easily adapted to concentra- 
tion of small elements in areas of interest and assignment of larger 
elements to other areas. 

The puirpose of this investigation is to determine the sensitivity of 
a FEM model to different grid geometries. The model used is a two- 
dimensional advection model with diffusion. It was deliberately chosen 
because of its simplicity. By keeping the model simple, all variations 
in the results can be attributed to differences in geometric and initial 
conditions rather than inaccuracies in the equation formulation. In 
addition to sensitivity studies, the model can be adapted as an opera- 
tional predictor of any scalar field such as cloudiness, vorticity, 
temperature or precipitation. 
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The Galerkin technique is employed to transform the differential 
equation for the dependent variable into a set of algebraic equations. 
Triangular elements are used with the area coordinate system. Later 
chapters will show that these techniques greatly simplify the 
computations . 
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II. THE METHOD 



The Galerkin technique will be explained by first an example/ after 
Crandall (1956) / and then the application to the FEM, after Kelly and 
Williams (1976) . 

A. EXAMPLE 

Consider a simple differential equation/ that which governs the cool- 
ing of an object 



with the solution defined on the interval 0 < t < 1. 

The critical step is to select a trial family of approximate solutions. 
(The members of a trial family are often called basis functions.) For 
simplicity/ consider the second order trial-solutions 



X = 1 



t = 0 




II-l 



t > 0 



2 

X = 1 + c_t + c^t 
1 2 



II-2 



/\ 



where x is the approximate solution 



The object is now to determine the coefficients c^ and c^. Next / form 
the residual/ R(t) from II-l 




II-3 
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* 



substituting II-2 into II— 3 



R 



1 + ^2 



(2t + t^) 



II-4 



Now we want to adjust and so that II-4 stays close to zero on the 
interval 0 < t < 1. 

The above discussion applies to any weighted residual technique; 
however^ Galerkin's technique has two requirements: 

1. The weighted coverages of the residual must vanish over the interval. 

2. The weighting functions must be the same functions of t as were used 

to construct the trial family (in II-2) . In this case those 

2 

weighting functions are t and t . 

That is, the weighted averages are formed by multiplying the residual 
by the weighting functions, integrating the product over the interval 
and setting the result equal to zero. 



1 




o 



t R dt = 0 



and 



/ 



t R dt = 0 



II-5 
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Substituting II-4 into II-5 



1 1 
t R dt = 



j t R dt = [t + 



(t+t^) + (2t^+t^) ]dt 



15 11 

= — + — c + — c =0 
2 6 1 12 2 



II-6a 



and 



j t^R dt = J [t^+ 



2 3 3 4 

(t +t )+ c^C2t +t )]dt 



o 

— + — c + — c =0 
3 12 1 10 2 



II-6b 



Solving II-6a and b for and we find: 
= -.9143, = 0.2857 



So that the approximate solution is 

X = 1 - 0.9143 t + 0.2857 



The exact solution to Equation II-l is 



-t 

X = e 
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Figure 1 shows the plots of the exact solution and the approximate solu- 
tion. Note that this technique yields the best fit only on the desired 
interval 0 < t < 1 and that the solutions may diverge from the interval. 

In general, Galerkin*s technique will reduce a partial differential 
equation to a family of N algebraic equations, where N is the mamber of 

basis functions. In this example Galerkin’s technique has reduced II-l 

2 

to two equations, II-6a and II-6b, in the two basis functions t and t . 

B. APPLICATION TO THE FEM 

The Finite Element Method divides the domain into discrete segments 
called elements with points (called nodes) arranged about the perimeter of 
the elements. The basis functions are then defined locally. These basis 
functions are usually low order polynomials that must be piece-wise 
continuous. A one-dimensional example is Figure 2 wherein the domain 
(x-axis) is divided into four elements (line segments) W through Z. All 
of the basis functions (A-C) are linear rising to a value of unity over 
each node (points 2 through 4) and are zero elsewhere. Now, if a variable 
S is defined over node 3, for example, then it may be approximated by 



S 



^2 A ^3 B U C 



II-7 



where Y^ Y^ values of S at nodes 2, 3 and 4 respectively 

and V^, V^, are the basis functions A, B, and C respectively. 

Equation II-7 may be substituted into any equation in which S appears. 
An equation using II-7 can be written for each node, by multiplying the 
basic equation by the basis function and integrating over the domain. 
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Figure 1, Solutions to cooling object example. Approximate solution 
is the curve connecting the *x* signs. Exact solution is 
the curve connecting the *+* signs. 
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Figure 2* One— dimensional example of linear basis functions. 





Figure 3. A domain with eight nodes and eight triangular elements. 
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The end result is a system of N equations in N unknowns, where N is the 
number of nodes. Fortunately^ in matrix form, these equations are often 
tri diagonal and may be solved with great economy. 

In two dimensions the basis functions are somewhat different, but the 
mathematics is identical. In Figure 3, the domain has been divided into 
eight triangular elements with eight nodes. Figure 4 shows the basis 
function (outlined in heavy black) at node 7 that must span all the 
elements about node 7. Note that the function has value of unity at 
node 7 and decreases linearly to zero at all surrounding nodes. Figure 5 
shows another basis function, but this one at node 5. Each node has a 
similar basis function about it. 

The value of a variable S may be approximated at any node i by 



where j ranges over all the nodes connected to node i including i itself. 



of the basis function of the i-th node. 

As in the one-dimensional case, an equation like II-3 is written 
for each node, i, then multiplied by the basis function and integrated 
over the domain. The resultant equation set, in matrix form, contains 
an n by n square matrix where n is the number of nodes. Although this 
matrix is usually not tridiagonal, careful selection of a nodal numbering 
system would yield a matrix that is strongly banded and fairly easily 
solvable. Alternatively, with these sparse matricies, iterative methods 
converge rapidly. 




j j 



V 



II-8 
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Figure 4. The linear basis function associated with node seven. 




Figure 5. Another linear basis function on the same domain. 
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Ill* EQUATION FORMULATION 



The equation of interest to the investigation is the two-dimensional 
advection equation with second order diffusion: 



3S 

3t 



+ 



3s 



3s 



3x ^ 3y 






v^s = 



0 



iii-i 



Define an inner product of two functions f(x,y) and gCx^y) as 




fg dx dy = <ffg> 



y X 



Following the Galerkin technique, form the residue and find the weighted 
averages, using the basis functions as weighting functions. Neglect the 
diffusion term at this stage. 




y X 

or 




V> + <u 



3s 

3x 



v> 



+ 




v> 



0 



III-2 



where V is the weighting function. Although the domain of integration 
is over the entire range of x and y, it will soon be shown that V is zero 
everywhere except locally around a specific node. So that for any node, 
III-2 is non-zero only when within one grid increment of that node. 
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The variables will have the form 



S = Y- V. 
3 3 



u = a. V. III-3 

3 3 



V = 6 . V. 
3 3 



where the repeated subscripts indicate sxarranation over the range of the 
sxibscript. The coefficient is a function of time only and the basis 
function is a function of space only. That is 



S = Yj Vj = 2 Yj (t) (x,y) 



etc. 

The critical step is the requirement that the weighting function in III-2 
be the same as the basis function in III-3. Substitution of III-3 into 
III-2 yields 



<Y .V. 
3 3 



3V. 

'V " "“kVj 7T 



v.> 

1 



+ 



<B V Y. 
k^j 




V.> 

1 



0 



III-4 



where the dot indicates a time derivative. 

Applying first the Galerkin technique and then the Gauss Divergence 
Theorem to the diffusion term yields 
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yy*K^(V^S)V^ dxdy ^ \ -{VS)V^ 

y X y X 



dxdy 









•VS) - vs • Vv^] dxdy III-5 



y X 



= K^{ V^VS • n dr - 



//"*' 

y X 



Vs dxdy} 



where n is a unit vector normal to the domain and dr is the differential 
distance along the path of integration on the perimeter of the domain, 

3s 

The contour integral in II 1-5 vanishes because cancels out due to 

3s 

cyclic continuity and = 0 because there is no flux across the channel 
walls at the north and south edges of the domain. Using III-3^ the re- 
maining term on the r.h.s. of III-5 can be written 



3v. 3v. 3V. 3v. 

-K^<W^.VS> = - , y. , y. 



III-6 



Combine III-4 and III-6 into the transformed advection equation as 
follows 



. 3v . 3v . 

Y.<v.,v.> + Y.-fa <v, - 5 -J- , v.> + s, <V, - 5 -^ , v.> 

3 3 ^ ]kkdx i k kdy i 

3v, 9V. 3v. 9v. 
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Now using a centered difference in time gives the forecast equation 



3V . 9v . 



3v. 3v. 3v. 3v. 



III-7 



when n is the time step* In this model and are not forecast 
quantities, but are either specified for each time step or are constant. 
They represent the prescribed wind field which is advecting the scalar S. 

This equation is valid for each node i and as such is a matrix 
equation of the form 

[A] {x} = {b} 



where 



[A] = < ^ 



with dimensions n by n 



U) 



< - v;- 



with dimensions n by 1 



{b} = The right-hand side of III-7 

with dimensions n by 1 

and n is the number of nodes. 

Note that all of the inner products are independent of time so that 
they need only be calculated once and placed in mass storage to be 
accessed at each time step. 
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IV. EVALUATION OF INNER PRODUCTS 



Even though the choice of elements in this investigation is triangles 

with nodes at the vertices^ this is by no means the only choice possible. 

The elements could be triangles with six nodes each; rectangles with six, 

eight or twelve nodes; or any other polygon or curved-sided element. More 

complex elements would allow higher orders of continuity. Bathe and 

Wilson (1976) describe the technique of isoparametric interpolation as 

applied to the FEM. However, these more complex elements require even 

more complex formulations. The mathematics can get out of hand very 

quickly requiring many more terms and complicated integrations. Also 

Galerkin's method is not the only weighted residual technique available. 

Crandall (1956) describes the collocation, subdomain, and least squares 

methods as well. Lower order continuity is deemed to be acceptable in 

this investigation because the grid resolution is fine enough so that 

the change in variable is assumed to be linear between grid points. 

This assumption allows the selection of triangular elements and the 

Galerkin method, each of which will significantly simplify the mathematics. 

Figure 6 shows a typical element with a point p described in area 

coordinates. Lines are drawn from p to each of the vertices (x . , y., 

3 3 

j=l,2,3) dividing the element into the areas A^ (j=l,2,3). 

Define: L. = A. /A 

3 3 

where A = total element area 
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Figure 6. A typical element illustrating the area coordinate 
system. 




Figure 7. An element with part of a basis function. 
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u. 



Now 



or, in matrix form 






X = l-i*! + ^ 2=^2 * ^’‘3 



y = + LjYj ■^ L3Y3 








or, solving for the L vector 



where 




2A 



" ^2 



^2 “ ^3 



^3 = ^2 - 



2 A 






2 A 


>=2 


^2 


2 A 


^3 


^3 


It 

r— 1 


^2 - 


^3 


It 

CN 


^3 - 


^1 


^3 = 


^1 “ 


^2 




IV- 1 



So that any one of the equations in IV-1 is really 



b . a . 

L . = 1 + X + — r y 

D 2A 2A 



IV-2 
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But by the chain rule 



and 



3x 



3l 

3x 



i 3 ^ ^ 3 

9l. 2A 3l. 

: J 




3y 3y 3 l^ 




2A 3l. 
3 



IV- 3 



where the repeated subscript indicates a suirmiation from one to three. 

Consider Figure 7. The heavy black lines outline an element 

divided into area coordinates defining point P, and only is labeled. 

The hatched triangle (labeled V^) is that part of the basis function 

associated with node 2 that lies over the element. There are two more 

basis functions called and that are associated with nodes 1 and 

3 respectively. Both of these functions also have linear sections over 

the element. However, they have been omitted from Figure 7 to avoid 

clutter. The altitude h^ is defined as unity. The altitude h is the 

value of at P. Recall ^ moves from node 2 to the 

opposite side of the element, h varies linearly from 1 to 0. Following 

the same path also varies linearly from 1 to 0. So that at any 

point P on the element L = h = V . In general, L. = V.. This is the 

z z j j 

first great advantage of triangular elements. 

So now, using IV- 3 



3v. 




3v. 




3V. 


b. 


3v. 


1 _ 


1 




+ - Z . 


1 




1 


3x 


2A 


3L^ 


^ 2A 


3l^ 


^ 2A 


3 L 3 




b. 


3l. 




3l. 


b^ 


3l. 




1 


1 


2 


1 


3 


1 


— 


2A 


3l^ 


2A 


3L^ 


2A 


\ 






1 




2 




3 
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] 



but 



So 




if i = j 

if i ^ j 




and similarly 



9v. a. 
^ 

3y “ 2A 



IV-4 



Note that these derivatives are dependent upon the individual elements 
and are constant with time. Therefore they need only be calculated once 
and be stored as part of the inner products. 

The second great advantage of triangular elements and linear basis 
functions is the general area integration formula 




mini 

(m + n + 2)1 



2A 



IV- 5 



or, rewritten as an example of an inner product: 



<L 



1 



2 






2111 

51 



2A 



A 

30 



Now the inner products for equation III-7 may be rewritten with the 
aid of IV-4 and IV-5: 
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if i ^ j 



<V., V> = <L^, !,.> = 2A = ^ 



1 41 6 



<V, 



3V. 

1 

k 3x 



b. 


b. 




b. 


J. 


1 

tl 

A 

V 


A 


_ D 


2A 


k' 1 2A 


12 
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b. 








_ D 
12 



9v. 

D 



a , 



''i" ■ 54 



if k 7 ^ i 



3v. 


3v. 


b. 


b. 




-J-> = 


1 




9x ' 


9x 


2A 


2A 


3v. 


9V. 


a. 


a . 


. 1 


—L> 


1 


-J. 


'9y ' 


9y 


2A 


2A 



if i = j 



if k 7 ^ i 



if k = i 
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V. INITIAL AND BOUNDARY CONDITIONS 



Simple initial conditions are used for the scalar field S and the 
wind components u and v. The S field consists of a half sine wave in 
the y direction and multiple waves in the x direction. These multiple 
waves are the mechanism used to generate smaller scale features. The 
term in the u equation causes the area of fastest flow to be coinci- 
dent with the area of highest spatial resolution. 



S = A sin (r 



2n7Tx 



L 




V-1 




V-2 



V = V = 0 



where 



A = arbitrary amplitude = 400 meters 



n = wave numbers 



L = channel length 



2400 km 



W = channel width 



2400 km 



U = mean zonal wind 



20 m/sec 



B = U X R; R is the ratio of the perturbation 
wind to the mean zonal wind 



Figures 8 and 9 depict the initial fields of S and u respectively. 




Figure 8. Initial field for the scalar S. 



o 




Figure 9. Initial u field with wind ratio R = 0.4. 
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Boundary conditions are also handled simply. Cyclic continuity is 
assured in the x-direction by selecting the proper node convectivity . 
Two techniques are used to control flow across the channel walls at 
the north and south sides. First, the exponent on the y-component of 
the initial wave is high enough to keep the value of S close to zero 
in the vicinity of the walls. Second, the solution method is iterative 
so that the values of the nodes one row in from the channel walls are 
set equal to the nodes on the walls at each iteration. The result is 
that diffusion from the edges into the interior is disallowed. 
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VI. COMPUTATIONAL TECHNIQUES 



A. NUMBERING SCHEMES 

As described in Chapter IV, the inner products may be easily 
evaluated on the individual element level. To facilitate this evalua- 
tion, a local numbering system is required for each element. An array 
called lEL, dimensioned N by 3, contains the identification number of 
the nodes on the three verticies for each of N elements. For each 
element these node numbers must be stored sequentially in a positive 
sense, that is counterclockwise. The node with which to start number- 
ing, however, is arbitrary. 

In addition to the local numbering system, a "mass” matrix is 
needed. The mass matrix is a matrix of coefficients whose rows are 
the equations of the system to be solved. There is one equation (row) 
for each node. Each equation will have a term (column) for each node. 

A non-zero entry in the mass matrix indicates connectivity. That is, 
if entry (I,J) is non-zero then node I is connected to node J. Two 
nodes are connected if they are both vertices of the same triangular 
element. In this model, each node is connected to a minimum of two and 
a maximum of six or eight other nodes depending upon model version. 

The model uses an array called MAME dimensioned N by 7 (or 9) that con- 
tains the numbers of all nodes connected to any specific node, includ- 
ing itself, for each of N nodes. I4AME will contain from 3 to 7 (or 9) 
non-zero entries in each row. 
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To save storage space, MAME is next compacted into a one dimensional 



vector NAME which contains only the non-zero entires of MAME. The 
utility vectors ISTART(N) and NUM(N) are also constructed in order to 
decode NAME. ISTART(N) contains the index in NAME of the starting 
position of the continuity of any node N. NUM(N) contains the length 
of the connectivity of node N within NAME. As an example of how to re- 
trieve the connectivity of a node, say node 14, from NAME; Let J = 
ISTART(14) and let K = NUM(14) . Then NAME(J) through NAME (J+K-1) contain 
the connectivity for node 14. 

Consider the matrix equation [A] {x} = {b} from Chapter III. It 
should be clear that [A] is the mass matrix mentioned above. It is a 
very sparse matrix containing the inner products <V^ , V^> for row j and 
column i. In this model [A] is compacted into the vector {h} using NAT4E 
as a lookup table. The result is that there is no longer a real matrix 
equation, but rather just the product of vectors 



{h} {x} = {b} 



VI-1 



B. SOLUTION TECHNIQUE 

An advantage of economy of space has been gained, but VI-1 now re- 
quires an iterative solution. The solution chosen is a standard Gauss- 
Seidel technique which requires a reversal of direction of iteration 
every pass. An average of 16 passes are needed for each time step. The 
solution is considered to have converged to its final value when: 
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X . 
1 



£-1 



< 10 



-6 



£ 



X. 

1 



for every node i 



with; £ = iteration number 



X. = solution for node i in x 
1 



More efficient solution schemes are available, but this one was 
selected because of its simplicity. Minimizing computation times is 
not of high priority to this investigation for two reasons: first, the 
recently installed IBM 3033 is a fast machine, and secondly, since an 
advection model doesn't contain gravity waves, stability can be main^ 
tained with very large time steps. As an example, a twelve hour 
forecast with a grid size of 200 km and a time step of one hour requires 
less than eight seconds of CPU time. Even the largest version of the 
model with longer forecast length, higher resolution and linkage to 
the plotting software takes less than eight minutes. 

C. ROBERT FILTER 

Haltiner and Williams (1980) discuss the advantages of time averag- 
ing. Averaging operators act as low-pass filters that tend to remove 
high frequency waves while having little effect on long-period waves. 

Even very weak filters will remove high frequency noise if used at every 
time step. In this model, the choice of filters is that of Robert (1966) . 
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First, assume that an average field already exists at time 

level n-1 well as the unaveraged field at time level n. Then 
the model uses its predictor equation to deteinriine the solution vector 
{x}« For each node, an unaveraged predicted value is computed as 



S 



n+1 



n-1 



+ X 



Next, the corrected values at time n are determined from 



S 

n 



n 



+ Y (S 



n+1 



2S + S J 
n n-1 



VI-2 



where Y is the averaging coefficient. Finally, is stored in place 

of S , and the model continues on to the next step. 

n-1 

The value of y must be carefully chosen as it does effect the com- 
putational stability. As y increases, the maximum time step decreases. 
Haltiner and Williams prefer a y less than 0.25 in order to permit a 
reasonable time step. Additionally, it should be noted that a large 
value of Y will begin to dampen waves in the meteorological frequency 
range. On the other hand, Robert used a value of y = 0,01 applied over 
a total forecast length of many days. This very weak filter was all 
that was needed for noise suppression. 

For the model in this investigation Y is an input parameter. For 
one set of tests in Chapter X, the value of y is varied from 0.01 to 
0.04. 
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D. DETERMINATION OF TIME STEP 



Cullen (1973) calculated the stability criterion for the two-dimen- 
sional problem as: 



At < 



Ax 

cv^ 



VI-3 



the minimum grid distance, Ax, varies from 200 km to about 40 km. The 
propagation speed, c, varies from 40 m/sec to 0. Accordingly, the time 
step has been chosen as 60 minutes for the coarse mesh versions, and 
ten or five minutes in the finer mesh versions. 
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VII. NODAL GEOI4ETRIES 



Once the decision has been made to use triangular elements, the next 
choice is the type of triangle. This investigation concerns both right 
and equilateral triangles. 

A. RIGHT TRIANGLES 

Figure 10 is an example of the right triangle arrangement used in 
this model for a uniform grid. The domain is divided into a series of 
rectangles, and then each rectangle is bisected by a diagonal to form a 
pair of triangles. Kelly and Williams found that significant bias was 
introduced into the model if all of the diagonals sloped in the same 
direction. This model overcomes that kind of bias by alternating the 
slope of the diagonals from rectangle to rectangle both horizontally and 
vertically. Consequently, the connectivity of the nodes will vary from 
four to nine total connections, and each node connects with itself. 

Two methods are used to vary the grid resolutions with right 
triangles. In the first method a FORTRAN DATA statement specifies co- 
efficients used in the calculation of the nodal coordinates. This method 
is very simple and has the advantage of keeping all the nodes in a 
rectangular arrangement. Rectangularly arranged nodes lend themselves 
well to harmonic analysis with some interpolation. The grid may be 
varied in either one or both the x and y-directions. Figure 11 is an 
example of the nodal arrangement by this method with variation in both 
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Figure 10. The domain divided into right triangular elements. Only 
the elements in the lower-left corner are drawn, but the 
pattern is repetitive over the entire domain. Each node 
is marked with an *X*. 
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Figure 11, Same as figure 10 with non-uniform elements. 
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directions. The disadvantage here is that the resolution must be kept 
relatively fine in parts of the domain away from the area of interest 
(the center) . Optimumly, the resolution should be fine about the area 
of interest and more coarse elsewhere, hence the next method. 

In the second method, the domain is divided into quadrants with the 
origin at the center of the domain. The nodes around the outer boundary 
of the first quadrant are specified as having the same coordinates as 
the unifom grid. But the nodes along the abscissa and ordinate are 
compressed toward the origin by the use of another DATA statement. Con- 
sider four nodes each with coordinates (x^, f i = 1...4 , Node 1 is 
on the ordinate, node 2 on the right boundary, node 3 on the abscissa, 
and node 4 on the top boundary. The line from node 1 to 2 is defined 
by 




And the line connecting nodes 3 and node 4 is 




X - X3 ^4-^3 



VII-1 



VII-2 



Equations VII-1 and VII-2 are then solved simultaneously to calculate the 
coordinates of the node located at the intersection of the two lines. By 
selecting the proper pairs of nodes along the boundaries and axes, the 
coordinates of all the nodes on the interior of the first quadrant are 
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found by VII-1 and VII-2, The nodes about the origin are then slightly 
adjusted to allow local uniformity* Finally the first quadrant is re- 
flected about the axes to form the other three quadrants in the domain. 
Figure 12 is an example of such an arrangement. 

The advantage of this method is that the resolution is fine only in 
the area of interest, and changes smoothly in all directions away from 
that area. The disadvantages are the complexity and the non-rectangu- 
larity of the nodes. 

B. EQUILATERAL TRIANGLES 

Preliminary results of the model indicated the presence of a degree 
of noise when using right triangles. Cullen (personal conversation) 
suggested the use of equilateral triangles as an alternative element 
arrangement. Figure 13 is an example of such an arrangement. The 
right triangles along the sides of the domain are in reality halves of 
equilateral triangles. Because of cyclic continuity, the domain actually 
"wraps around" so that all of the elements are equilateral. Note that 
the maximum connectivity is now only seven and that the domain is no 
longer square but rectangular. 

Variation of resolution is now somewhat more complicated and is 
accomplished by a transformation of coordinates. The desired end is a 
smooth and easily controllable stretchage and shrinkage of the coordi- 
nate axes. In the x-direction, the transformation is 

X = X + a cos kx VII-3 
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Figure 12. An alternate method of varying resolution. 




Figure 13. Same as figure 10 except with equilateral triangular 
elements. 
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where 



X = transformed coordinate 



X = original coordinate 
a = coefficient to be determined 




L = channel length 

9x 

the map factor is defined by 



tt-- = 1 - ak sin kx 

dX 



whose maximum and minimum values are 



(—) = 1 
^9x max 



+ ak 



(— ) 

3x min 



= 1 - ak 



The ratio, r^, of maximum stretch to minimum shrink is 

1 + ak 
^1 1 - ak 



or, solving for a 



a 




- 1 
+ 1 ) 



VII-4 



To keep the transformed coordinate system from folding back on itself, 
a must be positive. That is, r^, must be equal to or greater than one. 



Variation of resolution is accomplished by selecting an appropriate 
value for r^ and then substituting equation VII-4 into VII- 3. Similar 
expressions can be derived for the transformation of the y-coordinate 
by the use of the ratio r^. 

The advantage of this method is the extremely sensitive control upon 
resolution afforded by the ratios r^ and r^- Once programmed, this method 
is much easier to use than the more cumbersome DATA statements of the 
previous methods. The major disadvantage is that the transformed system 
is no longer made up of equilateral triangles if one of the ratios is 
different than one. Figure 14 is the nodal arrangement with - and 

^2 = '• 



0 



46 



I 



i 






2078 




Figure 14. Same as figure 13 but r^ - 5, r^ - 3* 
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Figure 15. Initial S field for the diffusion tests. 
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VIII • FORCING TERM 



As part of an effort to improve the testing of a numerical scheme, 
it is often beneficial to force the scheme with the desired solution, 
if known. Since this model is a simple advection scheme, the output will 
look very similar to the initial condition. 

Assume the solution is of the form 



S = cos [y(X - At) + 6] 



where 




A = phase speed. 



X = transformation of x coordinate 



If V = 0, then the forced advection equation is 




VIII-1 



but 



^ = yX sin[y(X - Xt) - j 



VlII-2 



and 




VIll-3 
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Substituting VIII-2 and VIII-3 into VIII-1 



y 7T 

y(X - u ~) sin[u(X - Xt) - j ] = F(x,t) 

which is equivalent to 

TT TT 

p(A - u -^) [sin(yx - — ) cos(yXt) - cos (yX - — ) sin(Xyt)] = F(x,t) 
9x 2 2 



applying Galerkin's technique this equation eventually reduces to: 



(f, - g.) <V., V.> = F(x,t) VIII-4 

: D D 1 

where 

9X TT 

f . = y(A - u TT') sin(yx - •-) cos yXt 

j 9x 2 ^ 

VIII-5 

g. = y (A - u — ) cos(yX - -r) sin yAt 

1 dx 2 



Now, VIII-4 should be included on the r.h.s. of the predictor equation 
III-7. Notice that both of the equations VIII-5 contain a time depen- 
dent part. The technique now is to calculate the time independent 
coefficients of VIII-5 at the onset and store them for the duration of 
program execution. During the forecast sequence, at each time step, 
the time dependent parts of VIII-5 are calculated, multiplied by the 
applicable pre-stored coefficients and then assembled into the mass 
matrix as indicated by VIII-4. 
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This forcing term will not be used in most of the model runs as a 
standard part of the predictor equation, but will be the subject of a 
specific set of tests designed to determine its contribution to 
forecast accuracy. 
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IX. ANALYTIC SOLUTION 



Consider the one dimensional equivalent of equation III-l with no 
diffusion: 



3t 



+ u(x) 



3x 



0 



IX-1 



where flow is zonal only and the velocity may vary spatially. Recall 
from Chapter V that the zonal wind equation is 



u(x) = U + B cos (kx - y) IX-2 

3tt 

where k = — / L = zonal wave length. 

L 

In a simple one dimensional advection equation for any time t, the 
following holds: 



S(x,t) = 

where x^ is some point upstream in the initial field. Now if x^ can be 
determined as a function of the current location and the elapsed time, 
then the analytic solution for the forecast at that time would be 
achieved. That is. 



Xq = F(x,t) 
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C - 

Since u = — 



then from IX-2 



^0 — 

-r— = U + B cos (kx^ 
at U 




Integrating 



f — %■/ 

J U + B cos (kx^ - j) J 



kdt 



With the help of a Table of Integrals, this expression is evaluated as 



(U-B)tan -^{kx - 
ttan [ — ] 



(U^ - 






, (U-B)tan i(kx- - |) 

- """ ^ -2 21/2 ' 

(U^ - B ) '^ 



Let 



r 



U 




-- B 



So that 



tan ^ [r tan j(kx - ] -tan ^ [r tan j(kx^ ~ j) ] = ~ kt 
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Solving for 

X* = ~ tan ^ {— tan [tan ^ (r tan --(kx - “ B^) kt] } IX-3 

0 2k k r 2 2 2 

The analytic solution at time t for any node at coordinates (x,y) 
is found by substituting x and t into IX-3 to obtain x^. Then x^ and 
y are substituted into equation V-1. This process requires each value 
of X to be operated upon by five trigonometric functions, each of which 
further compounds truncation error. As a result, this subroutine must 
be run in double precision. Even so, there is still some small error 
in the analytic solution field. This investigation is concerned with 
the relative errors of various nodal arrangements and is not intended 
for comparison studies between finite element and other methods. Since 
the errors in all test cases are relative to the same solution, the 
small error in the solution is not considered significant. 
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X, RESULTS 



The model contains eight variables: element type (equilateral and 

right triangles ) , east-west resolution^ north-south resolution, wave 
number, diffusion, Robert filter, uniformity of flow, and number of nodes 
on the domain. Additionally, the time step may be varied and the 
initial conditions may be changed so that sharper waves are created. 

If all of these variables were tested throughout their meteorological 
range, and if all the interactions between the variables were investi- 
gated, the number of computer runs required would be well over ten 
thousand. Since that many runs is not practical, decisions have to be 
made concerning the scope of this investigation. There seems to be one 
of two ways to proceed: 

1. One or two questions could be answered rather exactly with an 
attempt to quantify the relationship between only a couple of variables. 

2. General answers can be attempted for many more questions with 
the aim to set qualitative guidelines for future finite element studies. 

In the first alternative above, quantitative relationships would be 
of little value unless they could be extended to FEM models as a class. 
Proof that such an extension is valid is clearly beyond this investiga- 
tion. As a result, the second alternative is chosen and only rather 
general questions are qualitatively answered. 

The model is run a total of 137 times. The runs are compared and 
contrasted in a number of test cases. Several test cases comprise the 
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investigation of each question. Two types of statistics are calculated 
for each run. First of all, a harmonic analysis is accomplished using 
nodes along each latitude circle, for the initial condition, the forecast 
field and the analytic solution field. Since the analysis requires 
equally spaced values, there is an interpolator that calculates values 
at specific points. The linearity of the basis functions lends itself 
very well to linear interpolation. One subroutine determines in which 
element any particular interpolation point lies. Another routine uses 
the coordinates of that point and the values at the nodes of the proper 
element to calculate the interpolated value. In the runs that include 
the forcing term (see Chapter VII) , the interpolater is disconnected and 
the raw nodal values are used. This harmonic analysis generated ampli- 
tude and phase information for all possible wave numbers. In the follow- 
ing paragraphs the term**phase speed” denotes the ratio of the phase 
shift of the forecast wave to the phase shift of the wave in the analytic 
solution. 

The second type of statistics generated considers the forecast and 
solution values as an ordered pair for each node. A root-mean-square- 
error (RMSE) , correlation and bias are calculated for this set of pairs. 
By-products of these statistics that are not referred to hereafter are 
the applicable means and standard deviations. 

The following sections are organized such that each one addresses 
an individual problem, question or technique. To this end, within each 
section, one or more of the eight variables listed above is varied from 
a standard or nominal model configuration. Changes in the generated 



55 



statistics are then attributed to the influence of the altered variables. 
The nominal model configuration follows: equilateral elements, uniform 

grid (r^ = r^ = If see Chapter VII), wave number one, no diffusion, 

Robert filter, y = 0.1, wind pertixrbation ratio R = 0.0 (see Chapter V), 
nximber of nodes = 13 X 12. The standard forecast length is 48 hours, 
during that period any feature, moving with the mean flow, should move 
3456 km or around the domain 1.44 times. Section A is a general statement 
based upon nearly all of the runs. It includes changes made to all of 
the variables. In all of the remaining sections, the departure from the 
nominal configuration is annotated. 

A. ACCURACY VS. RESOLUTION 

For this test, six runs are discarded because they are deliberately 
constructed to test unique factors discussed later. In the other 131 
runs, the forecast values stay within reasonable limits for all combina- 
tions of resolution. In general, the accuracy over the domain as a whole 
decreases as the stretch of the grid resolution increases, but the 
statistics all remain bounded. The correlations are all above 0.96, the 
phase speeds are within four percent and the RMSE is always less than 

eight percent of the range of the amplitude. Typical values are cor- 
relation 0.99, phase speed 1.008 and RMSE one percent of range. 

This result is very encouraging in that no cases have to be dis- 
regarded simply because they can not be explained. The cases are 
consistent with one another and exceptions are few. 
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B. EFFECT OF DIFFUSION 



A particularly noisy case is used for this test. The domain length 
scale has been doubled so that the non-uniform wind is more obvious. 

With resolution of r^ = 1, r^ = 2 wind perturbation ratio R = 0.4, and 
13 X 24 nodes, the diffusion values are incrementally increased from 
zero. A total of eight runs are made. Figure 15 is the initial condi- 
tion for this test and Figure 16 is the 48 hour forecast with no diffu- 
sion. As diffusion is increased the reduction in RMSE is dramatic and 
the correlation increases, but the wave amplitude diminishes and the 
maximum time step must be decreased to insure stability. However, a 
critical point is eventually reached at which the amplitude is reduced 
so much that the RMSE starts to climb again and correlation falls. 

Figure 17 is the 48 hour forecast with the value of diffusion that 
nearly minimizes the RMSE. Note the absence of noise but the obvious 
loss of amplitude. Compare this with Figure 18, the analytic solution. 
Notice that the contour packing is downstream of the high in Figures 16 
and 18, but upstream of the high in Figure 17. This reversal of packing 
is a by-product of the manner in which the model applies diffusion. 
Cullen (personal communication) suggests that diffusion be increased in 
the areas of fastest flow and decreased where the flow is slower. 
Following this advice, the model applies diffusion proportional to the 
windspeed at each node. The effect with non-uniform flow is that diffu- 
sion is also non-uniform, hence the slight alteration in the forecast 
field. Diffusion appears to be a very effective noise filter, but tests 
must be made independently for each configuration because the useful 
limit of diffusion will vary with nodal geometry. 
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Figure 16. 48 hour forecast for field in figure 15 with no diffusion. 




Figure 17. Same as figure 16 but with optimum diffusion. 



58 



2078 



6 

>» 



0 




Figure 18. Analytic solution for diffusion tests. 




Figure 19. Initial S field for Robert filter tests. 
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C. EFFECT OF ROBERT FILTER 



For this test the model is configured nominally to maximize the accu- 
racy of the solution. The only variable that is changed is the Robert 
averaging coefficient. Figure 19 is this initial configuration. The 
averaging coefficient, y starts at a value of 0.4 and is slowly decreased 
to 0.01 over 11 runs. The RMSE decreases from 8.7 to 1.74 over this range 
with a corresponding increase in correlation and the phase speed approaches 
unity. 

Since the filter is desirable to control high frequency noise, a com- 
promise must be reached that will not unduly degrade the forecast. For 
the rest of this investigation a coefficient of y = 0.1 is used since the 
associated RMSE is only 2.54. Figure 20 is the 48 hour forecast field 
with this coefficient and Figure 21 is the corresponding analytic solution. 
There is little improvement in RMSE left to be realized if any effective 
filtering is desired. 

D. ACCURACY VS. WAVE NUMBER 

Again the model is configured nominally, except that the number of 
nodes is 13 X 24. Six runs are made with the wave number vairying from 
one to eight. One would expect the forecast to deteriorate as the wave 
number increases. This is due in part to the fewer number of nodes along 
a latitude circle available to resolve each wave. But the deterioration 
of RMSE is primarily the result of the displacement error measured in 
terms of the wavelength. A small displacement error with a short wave 
would produce a much greater deviation from the true solution than would 
be produced by the same displacement error in a long wave. 
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20. 48 hour forecast for field in figure 19- 




Figure 21- Analytic solution for Robert filter tests. 
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The results of this test verifies this reasoning. As the wave number 



of the initial field increases the RMSE worsens, even though the phase 
speed error stays within four percent. Four of these cases were excluded 
from Section A above. 



E. RIGHT TRIANGLES VS. EQUILATERAL TRIANGLES 

Five pairs of runs are used to test the hypothesis that elements 
formed by equilateral triangles yield forecasts superior to those based 
upon elements formed by right triangles. In all cases the flow is 
uniform but the resolution is varied among the cases from uniformity 
to variations in either or both directions. Partial results appear in 
Table I as Cases 1 through 5. Figures 22, 23 and 24 are the initial 
field, 48 hour forecast and analytic solution, respectively for the 
right triangles in Case 3. Figures 25, 26 and 27 are the corresponding 
figures for equilateral triangles also in Case 3. 

The differences between the two arrangements, while not dramatic, 
are significant and consistent. In every case the RMSE is lower and 
correlation higher for equilateral triangles. The average reduction in 
RMSE is about 20 percent. Not shown in Table I but readily apparent in 



Wave Number 



Phase Speed 
Ratio 



1 

2 

3 

4 

5 
8 



1.008 

1.008 

1.02 

1.03 

1.04 
1.03 
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Figure 22. Initial S field for Case 3 with right triangles. 



o 




Figure 23. 48 hour forecast for field in figure 22. 
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Figure 24, Analytic solution corresponding to figure 23, 




Figure 25. Initial S field for Case 3 with equilateral triangles. 
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Figure 26. 48 hour forecast for field in figure 25. 




Figure 27. Analytic solution corresponding to figure 26. 
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the harmonic analysis, is the generation of a spurious wave number five 



in the right triangle arrangement in four of the five cases. 



TABLE I 

RESULTS OF RIGHT VS. EQUILATERAL ELEMENT 
AND DIRECTION OF FLOW TESTS 









RIGHT 


EQUILATERAL 


CASE 




^2 


RMSE 


CORR. 


RMSE 


CORR. 


1 


1 


1 


4.25 


.9994 


2.54 


.9999 


2 


2 


1 


4.54 


.9993 


4.17 


.9996 


3 


1 


2 


6.40 


.9991 


4.87 


.9996 


4 


2 


2 


6.96 


.9990 


5.60 


.9994 


5 


3.55 


1 


6.19 


.9981 


5.75 


.9992 


6 


1 


4 


7.48 


.9980 






7 


4 


4 


6 . 60 


.9990 






8 


1 


3.55 


9.00 


.9980 







F. EFFECT OF ADDING NODES 

Next, 16 runs are combined to form nine test cases in order to 
determine the effect that adding nodes has on the accuracy of the fore- 
cast. In this section, the variables are wave number, resolution and 
number of nodes. It is found that where the model does very well (long 
waves, uniform flow, relatively small difference in resolution over the 
entire domain) , the addition of more nodes does not increase the accuracy, 
especially if the additional nodes destroy the original equilateral 
triangularity. However, if the nodes are added symmetrically in both 
the north-south and the east-west direction, so that the effect is to 
increase the resolution uniformly while preserving equilateral triangu- 
larity, then the result is an improved forecast. 
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On the other hand, when the model is initialized with shorter waves. 



more nodes will improve the forecast even though they may not be added 
symmetrically. The major advantage of additional nodes is that they 
allow the resolution to be changed more slowly and more smoothly while 
attaining a desired minimum resolution. The arrangement with fewer nodes 
requires a faster and more abrupt resolution change in order to attain 
the same minimum resolution. The smoothness of the transition between 
maximum and minimum resolution is apparently the critical factor. Later 
paragraphs also address this factor. The phase speed ratios for the 
uniform grids are as follows: 



Nodes 


Wave 

Number 


Phase Speed 
Ratio 


RliSE 


13 X 12 


1 


1.001 


2.5 


13 X 24 


1 


1.008 


10.4 


25 X 24- 


1 


1.003 


7.4 


13 X 12 


2 


1.003 


60.29 


13 X 24 


2 


1.002 


30.7 


25 X 24 


2 


1.003 


25.27 



G. VARIATION OF RESOLUTION VS. DIRECTION OF FLOW 

The initial conditions for the model are deliberately set to allow 
only zonal flow. (See Chapter V.) One of the reasons for this is to 
be able to specifically test the effects of variation of resolution along 
the direction of flow. A total of 13 runs form six test cases. Many of 
these runs are the same ones as used in tests in preceding paragraphs. 
Some of these cases appear in Table I. In that table the comparisons 
are now vertical. For example, compare Cases 2 with Case 3 to determine 
that r^ = 2 and r^ = 1 is superior to r^ = 1 and r^ = 2 for both right 
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and equilateral triangles. This means that varying the resolution with 
the flow is superior to varying across it. Case 2 is illustrated in 
Figures 28 and 29 which are the initial field and 48 hour forecast re- 
spectively. They should be compared with Figures 25 and 26 which, as 
stated previously, are part of Case 3. The same result holds for every 
case on the table. Additional tests have also been run with resolution 
ratios of six and eight. The results are the same. 

One further test can be conducted with the data on Table I. Resolu- 
tion ratios of r = r = 2 will yield the same overall resolution as one 

J. ^ 

of the ratios equal to 3.55 and the other unity. This is the rational 
for Cases 4, 5 and 8. These cases show that better results can be ob- 
tained by varying resolution in both directions as opposed to only across 
the flow by an equivalent amount. This fact agrees with the results of 
the previous section in that varying resolution in both directions tends 
to maintain the equilateral triangularity better than varying across the 
flow alone. 

H. SMOOTH VS. ABRUPT VARIATION OF RESOLUTION 

For this test, a version of the model is employed that allows the 
resolution to be changed abruptly. That version is paired with a more 
standard version that allows smoother variation of resolution. Both 
versions use right triangular elements. Both versions generate the same 
maximum resolution, but with very different transitions to the area of 
coarser mesh. Twelve runs are paired into six tests. Two different 
wave numbers and three kinds of flow are tested. These runs included 
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Figure 28. Initial S field for Case 2. 




Figure 29. 48 hour forecast for Case 2. 
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both the normal and '* flipped" (See Section J below) configuration. The 
results are contained in Table II. Even though only six tests are used, 
the results are dramatic and highly significant. The reduction in RMSE 
from the abrupt to the smooth cases averages 62 percent. Although not 
shown, the reduction in bias averages 70 percent. Additionally the 
abrupt cases generate significantly more noise at all frequencies. Cases 
10 and 12 are the last two cases excluded from Section A above. Figures 
30 through 33 illustrate the evolution of the forecast in Case 9. The 
figures are the initial field and 12, 24 and 48 hour forecast respective- 
ly for the abrupt change. Compare these to Figures 34 and 35 which are 
the initial field and 48 hour forecast for the corresponding smooth change. 

Combining this result with that of Section F above, it is obvious 
that the model is highly sensitive to the rate and smoothness of the 
transition from fine to coarse resolution. 



TABLE II 

RESULTS OF SMOOTH VS. ABRUPT TESTS 





WAVE 


WIND 


SMOOTH 


ABRUPT 


CASE 


NUMBER 


RATIO 


RMSE 


CORR. 


RMSE 


CORR 


9 


1 


0.0 


4.54 


.999 


13.6 


.989 


10 


2 


0.0 


28.3 


.956 


89.2 


.470 


11 


1 


0.2 


9.46 


.996 


23.5 


.968 


12 


1 


0.4 


21.4 


.979 


39.9 


.895 


IIF 


1 


0.2 


4.53 


.999 


14.6 


.987 


12F 


1 


0.4 


8.83 


.996 


20.0 


.979 



I. EFFECT OF VARIABLE WINDS 

In several of the previous tests, reference has been made to non- 
uniform flow. For this section and the following two sections a total 
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Figure 30. Initial S field with abrupt change of resolution 
(Case 9) . 



o 




Figure 31. 12 hour forecast for field in figure 30. 
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Figure 32. 24 hour forecast for field in figure 30. 



o 




Figure 33. 48 hour forecast for field in figure 30. 
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Figure 34. Initial S field with smooth change of resolution 
(Case 9) . 



o 




Figure 35. 48 hour forecast for field in figure 34. 
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of 42 runs are made to study the effects of wind and the forcing term. 
The 42 runs are formed into 14 tests in which the wind ratios are 
specified as 0.0 for uniform flow and 0.2 and 0.4 for variable flow. 

The tests include cases that have both uniform and non-uniform resolu- 
tion, both equilateral and right triangular elements and the 13 x 12 
and 13 x 24 nodal arrangements. Table III contains the results of 12 
of the runs which comprise four tests. In 12 of the total of 14 tests, 
the model performs better with the uniform wind. In all tests once 
the wind is non-uniform the forecasts are degraded as non-uniformity 
increases. 



TABLE III 

RESULTS OF VARIABLE WIND, "FLIPPED" WIND AND FORCING TERM TESTS 



UNFORCED FORCED 



CASE 


RATIO 


RMSE 


SPEED 
. .NORMAL. , 


CORR. 


RMSE 


SPEED 


CORR. 


13 


o 

• 

o 


2.9 


1.004 


.9996 


2.9 


1.004 


.9996 


14 


CM 

• 

O 


3.4 


1.005 


.9995 


3.3 


1.005 


.9995 


15 


0.4 


CD 

• 

-J 


1.006 


.9970 


8.7 


1.006 


.9970 



FLIPPED 



16 


0.0 


00 

• 

CM 


1.005 


.9996 


2.8 


1.007 


.9996 


17 


0.2 


2.7 


1.004 


.9996 


2.6 


1.004 


.9997 


18 


0.4 


3.1 


1.005 


.9995 


2.9 


1.004 


.9960 



J. "FLIPPED" WINDS 

Up until now, the initial conditions, in conjunction with the trans- 
formation of coordinates scheme, require that the area of highest wind 
is also the area of highest resolution, provided, of course, that both 
the wind and the resolution are non-uniform. Now, downstream of the area 
of highest winds is an area of convergence and downstream of the area 
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of lowest wind is an area of divergence. Therefore, where the conver- 
gence runs into the divergence is the area where the gradients are 
highest. It is a "bottleneck". But that is exactly the area where the 
resolution is the most coarse. The forecast should improve if the grid 
or initial wind pattern could be "flipped" to concentrate the nodes in 
the area of minimum winds. Since "flipping" the grid requires fewer 
changes, that is the choice here. As an example. Figures 36 and 37 are 
the initial field and 48 hour forecast of such an arrangement. 

To verify this reasoning, 24 runs are formed into 12 tests. Cases 
13, 14, and 15 on Table III are "normal" whereas Cases 16, 17 and 18 
are "flipped". In all tests, the "flipped" arrangement yields superior 
forecasts. This result justifies the reasoning that the high resolution 
should be concentrated in area of highest gradients, not the areas of 
fastest flow. 

K, FORCING TERM 

The forcing term was developed in Chapter VIII as an alternate method 
of testing the forecast if the general form of the end product is known. 

The 42 runs form 21 pairs of forced/ unforced tests, some of which are in 
Table III. The forced forecasts are slightly better in 15 tests, the 
unforced are slightly better in five tests. One test comes out a tie. 

In general, the unforced forecasts are better if the wind is uniform while 
a non-uniform wind favors the forced forecast. However, all differences 
are very slight. Figures 38 and 39 are the initial field and 48 hour fore- 
cast respectively for a typical forced forecast with a uniform wind. They 
should be compared with Figures 28 and 29 which form the corresponding 
unforced case. 
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Figure 36. Initial S field with "flipped" winds, r 
and R = 0.2. 



2, r. 
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Figure 38. Initial S field for a typical forced case. 




Figure 39. 48 hour forecast for field in figure 38. 
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XI. CONCLUSIONS 



This report applies the Finite Element Method to the two-dimensional 
advection equation with diffusion. Although the underlying equation is 
simple enough, FEM models have certain inherent advantages, which include 
very accurate phase speeds and the ability to change resolution easily. 
However, they also have limitations, such as a relatively high computa- 
tional cost, and a requirement for fairly sophisticated and efficient 
programming. The purpose of this investigation is to establish some guide- 
lines and provide qualitative methods that may be used in the development 
of future FEM models. 

The conclusions are given below in the same order as they are discussed 
in Chapter X; 

1. The model yields reasonable and consistent results. There appear to 
be no major algorithmic or coding problems as all the results are 
mathematically and meteorologically plausible. 

2. Diffusion is an effective noise filter; however it does have a maximum 
usable limit which must be independently determined for each applica- 
tion. Numerical models should contain at least some provision for 
diffusion. 

3. The Robert filter is also a valuable feature to have built into a 
model. In general, the filter should be no heavier than absolutely 
needed to control high frequency noise. This model uses a filter of 
0.1; higher values are not recommended. 
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4. The model handles wave numbers 1 and 2 well (24 and 12 nodes per 
wave, respectively) and wave number 3 (8 nodes per wave) with less 
accuracy. Higher wave numbers maintain accurate phase speeds 
although RMSE deteriorates. 

5. Formulations with elements as equilateral triangles outperform those 
with elements as right triangles. The difference in RMSE is about 
20 percent. The right triangle formulation generates spurious short 
waves which may be be the primary cause of the relative inaccuracy. 

6. Adding nodes (and hence, elements) to the domain will improve the 
model if they are added symmetrically. More nodes allow smoother 
transition between fine and coarse resolution. This smoothness is 
critical. 

7. An unexpected result is that this model performs better when the reso 
lution is varied with the flow rather than across it. The reasons 
for this are unclear. 

8. The most important conclusion of this investigation is that the 
smoother and slower the change in the resolution, the better the 
forecast. In this model the smooth change reduces the RMSE by 62 
percent over the abrupt change. 

9. The model works best with uniform flow. In general, as the flow 
departs from uniformity, the forecast continues to degrade. 

10. The forecast can be improved by concentrating the high resolution 
areas in the region of the strongest gradients. 

11. Inclusion of a term that forces the supposed solution back into the 
model may yield slightly better forecasts. It forms a valid, alterna 
tive method to test numerical models. 
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